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We develop, test and compare new numerical and geometrical methods for improving the accuracy 
of extracting waveforms using characteristic evolution. The new numerical method involves use of 
circular boundaries to the stereographic grid patches which cover the spherical cross-sections of the 
outgoing null cones. We show how an angular version of numerical dissipation can be introduced 
into the characteristic code to damp the high frequency error arising form the irregular way the 
jy^ , circular patch boundary cuts through the grid. The new geometric method involves use of the 

Weyl tensor component ^4 to extract the waveform as opposed to the original approach via the 
Bondi news function. We develop the necessary analytic and computational formula to compute the 
0(l/r) radiative part of $4 in terms of a conformally compactified treatment of null infinity. These 
methods are compared and calibrated in test problems based upon linearized waves. 
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PACS numbers: 04.20Ex, 04.25Dm, 04.25Nx, 04.70Bv 



I. INTRODUCTION 



The unambiguous geometric description of gravitational radiation in curved spacetimes traces back to the work 
of Bondi [l| et al, Sachs [|J and Penrose 0]. By formulating asymptotic flatness in terms of characteristic hyper- 
surfaces extending to infinity, they were able to reconstruct, in a nonlinear geometric setting, the basic properties 
of gravitational waves which had been developed in linearized theory on a Minkowski background. The major new 
nonlinear features were the Bondi mass and news function, and the mass loss formula relating them. This approach 
has been implemented as a characteristic evolution code [3, HI which computes the radiation field at infinity by using a 
Penrose compactification of the space-time. The code computes the gravitational radiation reaching infinity in terms 
of boundary data supplied on an inner worldtube. This has timely application to the important astrophysical problem 
of the inspiral and merger of a binary black hole. Several Cauchy codes, using an artificial outer boundary condition, 
are now able to simulate this binary problem. By using the data on a worldtube carved out of these binary black 
hole spacetimes obtained by Cauchy evolution, the characteristic code can supply the resulting waveform at infinity. 
In this work, we develop and test new methods designed to enhance the accuracy of this approach to computing 
■ gravitational waveforms, which has been called Cauchy- characteristic extraction (CCE) Q. 

The Cauchy codes presently being applied to the binary black hole problem introduce an artificial outer boundary, 
where some boundary condition must be employed. The choice of the proper boundary condition for an isolated 
radiating system is a global problem, which can only be treated exactly by an extension of the solution to infinity, e.g. 
by conformal compactification. The most elegant such approach is the extension of the Cauchy problem to future null 
infinity 2 + by means of a hyperboloidal time foliation [3] (see [1, Q for reviews of progress in this direction) . Another 
approach is to extend the solution to T + by matching the interior Cauchy evolution to an exterior characteristic 
evolution, i.e. Cauchy-characteristic matching (CCM) figl. ( see for a review). CCM has been applied successfully 
to gravitational wave computations in the linear regime [12| but has not yet been extended to the nonlinear binary 
black hole problem. 

When an artificial finite outer boundary is introduced there are two broad sources of error: 

• The outer boundary condition 

• Waveform extraction at an inner worldtube. 



The first source of error stems from the outer boundary condition, which must lead to a well-posed constraint- 
preserving initial-boundary value problem. This has not yet been fully established for any of the present black 
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hole codes. But, even were such boundary conditions implemented, the correct boundary data must be prescribed. 
However, this boundary data can only be exactly determined, in general, by extending the solution to infinity [T3 |. 
Otherwise, the best that can be done is to impose a boundary condition for which homogeneous boundary data, i.e. 
zero boundary values, is a good approximation. One proposal of this type [l4| is a boundary condition that requires 
the Newman-Penrose (l5j Weyl tensor component to vanish. In the limit that the outer boundary goes to infinity 
this outer boundary condition becomes exact. In the present state of the art of black hole simulations, this approach 
comes closest to a satisfactory treatment of the outer boundary [l6j] . 

The second source of error arises from waveform extraction at an inner worldtube, which must be well inside the 
outer boundary in order to isolate it from errors introduced by the boundary condition. There the waveform is typically 
extracted by a perturbative scheme based upon the introduction of a background Schwarzschild spacetime. This has 
been carried out using the Regge-Wheeler-Zerilli [1?], [3 treatment of the perturbed metric, as reviewed in 
and also by calculating the Newman-Penrose Weyl component ^4, as first done for the binary black hole problem 
in [2(1 HH, HH . In this approach, errors arise from the finite size of the extraction worldtube, from nonlinearities 
and from gauge ambiguities involved in the arbitrary introduction of a background metric. The gauge ambiguities 
might seem less severe in the case of "J 4 (vs metric) extraction, but there are still delicate problems associated with 
the choices of a preferred null tetrad and preferred worldlines along which to measure the waveform (see [24] for an 
analysis). 

Cauchy-characteristic extraction, which is one of the pieces of the CCM strategy, offers a means to avoid this error 
introduced by extraction at a finite worldtube. In CCE, the inner worldtube data supplied by the Cauchy evolution is 
used as boundary data for a characteristic evolution to future null infinity, where the waveform can be unambiguously 
computed by geometric methods. By itself, CCE does not use the characteristic evolution to inject outer boundary 
data for the Cauchy evolution, which can be a source of instability in full CCM. Highly nonlinear tests in black 
hole spacetimes Q have shown that characteristic evolution is a stable procedure which provides the geometry in 
the neighborhood of null infinity up to numerical error; and tests in the perturbative regime (25j show that CCE 
compares favorably with Zerilli extraction and has advantages at small extraction radii. However, in astrophysically 
realistic cases which require high resolution, such as the inspiral of matter into a black hole [26| . this error has been 
a troublesome factor in the postprocessing of the numerical solution which is necessary to compute the asymptotic 
quantities determining the Bondi news function. 

There are two distinct ways, geometric and numerical, that the accuracy of this calculation of the gravitational 
waveform might be improved. In the geometrical category, one option is to compute ^4 instead of the news function 
as the primary description of the waveform. We discuss this in Sec. IIII1 where we develop the extensive formulae 
necessary to compute the asymptotic 0(1 /r) part of ^4, i.e. \&4 = lim^oo r^4, which governs the radiation. 

In the numerical category, some standard methods for improving accuracy, such as higher order finite difference 
approximations, would be easy to implement whereas others, such as adaptive mesh refinement, have only been tackled 
for ID characteristic codes [27]. But beyond these methods, a major source of error in characteristic evolution arises 
from the intergrid interpolations arising from the multiple patches necessary to coordinatize smoothly the spherical 
cross-sections of the outgoing null hypersurfaces. The development of grids smoothly covering the sphere has had a 
long history in computational meteorology that has led to two distinct approaches: (i) the stereographic approach in 
which the sphere is covered by two overlapping patches obtained by stereographic projection about the North and 
South poles [28|; and (ii) the cubed-sphere ap pro ach in which the sphere is covered by the 6 patches obtained by a 
projection of the faces of a circumscribed cube [291 ]. Recently, the cubed-sphere approach has received much attention 
because the simple structure of its shared boundaries allows a highly scalable algorithm for parallel architectures. 
A discussion of the advantages of each of these methods and a comparison of their performance in a standard fluid 
testbed are given in (28j . In numerical relativity, the stereographic method has been reinvented in the context of the 
characteristic evolution problem [30| ; and the cubed-sphere method reinvented in the context of building an apparent 



horizon finder [31( . The characteristic evolution code was first developed using two square stereographic patches, each 
overlapping the equator. We consider here a modification, based upon the approach advocated in [28|, which retains 
the original stereographic patch structure but shrinks the overlap region by masking a circular boundary near the 
equator. Recently, the cubed-sphere method has also been developed for application to characteristic evolution [32l.[33|. 
These geometric and numerical considerations lead to four options for improving CCE: 

• Computation of the news function using circular stereographic patches 

• Computation of the Weyl tensor using circular stereographic patches 

• Computation of the news function using the cubed-sphere patching 



• Computation of the Weyl tensor using cubed-sphere patching 
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We compare these options here in the context of model problems designed to test their application to CCE. 
Because the cubed-sphere approach requires further code development to be applied to CCE, in Sec. [V] we present a 
test based upon the propagation of a wave on the sphere to provide a preliminary comparison with the stereographic 
approach. The test compares their accuracy in calculating the angular derivatives required in the news and Weyl 
tensor extraction algorithms. In Sec. IVI1 we next present tests of CCE which compare the news and Weyl tensor 
approaches in a linearized gravitational wave test problem. 

The development of finite-difference evolution algorithms, which was largely motivated by application to computa- 
tional fluid dynamics (CFD). It has utilized the method of lines, where a 3-dimensional spatial domain is discretized 
to yield a set of coupled ordinary differential equations in time for the grid values, which are then integrated, e.g. 
by a Runge-Kutta procedure. This 3 + 1 approach is not applicable to the 2 + 1 + 1 format of characteristic evolu- 
tion considered here, in which the discretization of a 2-dimensional spherical set of characteristics leads to coupled 
2-dimensional partial differential equations in the plane spanned by the outgoing and ingoing characteristics. This 
2 + 1 + 1 approach is natural to general relativity since the characteristics (light rays) are fundamental to describing the 
dynamical geometry of space-time. It would be impractical in CFD in which the characteristics have a complicated 
dynamic relation (determined by equations of state) to the fixed Euclidean geometry. As a result, characteristic evo- 
lution algorithms were developed only recently in the context of general relativity and there has been relatively little 
analysis of their computational properties. In particular, for CFD or any symmetric hyperbolic system, numerical 
dissipation can be added in the standard Kreiss-Oliger form (34[. One of the main results of this paper is to show how 
analogous dissipation can be successfully applied in a 2 + 1 + 1 format. In the original version of the PITT code, which 
used square stereographic patches with boundaries aligned with the grid, numerical dissipation was only introduced 
in the radial direction [35[ ■ This was sufficient to establish numerical stability. In the new version of the code with 
circular stereographic patches, whose boundaries do not fit regularly on the stereographic grid, numerical dissipation 
is necessary to control the high frequency error introduced by the intergrid interpolations, as previously noted in the 
treatment of a fluid problem using circular stereographic patches [28j j . We begin with a brief review of the formalism 
underlying the characteristic evolution code in Sec. [IT] and show how the essential new feature of angular dissipation 
can be incorporated. 

The two spherical grid methods, stereographic and cubed sphere, are briefly described in Sec. HVl We present the 
test results in Sees. IVl and IVII and we summarize our conclusions in Sec. IVIII 



II. CHARACTERISTIC FORMALISM 



The characteristic formalism is based upon a family of outgoing null hypersurfaces, emanating from some inner 
worldtube, which extend to infinity where they foliate I + into spherical slices. We let u label these hypersurfaces, 
x A (A — 2,3) be angular coordinates which label the null rays and r be a surface area coordinate. In the resulting 
x a = {u,r,x A ) coordinates, the metric takes the Bondi-Sachs form [H, Q 

ds 2 = - { e 2fi — - r 2 h AB U A U B ) du 2 - 2e 2f3 dudr - 2r 2 h AB U B dudx A 



r 

+ r 2 h AB dx A dx B , (2.1) 

where h AB h B c — &c and det(h AB ) — det(q AB ), with q AB a unit sphere metric. In analyzing the Einstein equations, 
we also use the intermediate variable 

Q A = r 2 e- 2 ^h AB U B . (2.2) 

The code introduces an auxiliary unit sphere metric q AB , with associated complex dyad q A satisfying q AB = 
\ (qaQb + QaQb) ■ For a general Bondi-Sachs metric, h AB can then be represented by its dyad component J = 
h AB q A q B /2, with the spherically symmetric case characterized by J = 0. The full nonlinear h AB is uniquely deter- 
mined by J, since the determinant condition implies that the remaining dyad component K = h AB q A q B /2 satisfies 
1 = K 2 — J J. We also introduce spin- weighted fields U — U A q A and Q = Q A q A , as well as the (complex differential) 
operators 3 and 8 (36|. Refer to d, H0| for further details. 

In this formalism, the Einstein equations G M „ = decompose into hypersurface equations, evolution equations and 
conservation conditions on the inner worldtube. As described in more detail in |37l. |38| . the hypersurface equations 
take the form 



13. r = Nf), 

U, r = r- 2 e 2p Q + Nu, 



(2-3) 
(2.4) 
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where [3C 



(r 2 Q), r = -r 2 (BJ + dK), r + 2r 4 d(r- 2 l3) r + N Q , (2.5) 
V r = \e^Tl - e^Be + ~r" 2 (r 4 (W + BU)) r + N w , (2.6) 

n = 2K- mx + -m 2 j + g 2 J) + — (§jgj - §jsj) (2.7) 

2 \K 

is the curvature scalar of the 2-metric h ab ■ The evolution equation takes the form 

2(rJ) jUr -(r- 1 V(rJ) }r ) r = 

-r- 1 (r 2 W) r + 2r- 1 e l3 B 2 e p - (r^V) r J + Nj, (2.8) 

where, Np, Njj, Nq, N\y and Nj are nonlinear terms which vanish for spherical symmetry. Expressions for these 
terms as complex spin- weighted fields and a discussion of the conservation conditions are given in [g] . 

The characteristic evolution code implements this formalism as an explicit finite difference scheme. In this paper 
we use second order accurate finite differences and we reduce all angular derivatives to first order by the introduction 
of auxiliary variables, as described in [3^ | 

A. Angular dissipation 

It is a feature of the composite mesh technique that numerical dissipation is necessary to stabilize the error intro- 
duced by intergrid interpolations. In the case of a square stereographic patch, whose boundary aligns with the grid 
lines, the dissipation built into the characteristic radial integration scheme is sufficient for this purpose [35j. However, 
because a circular boundary fits into a stereographic grid in an irregular way, angular dissipation is also necessary in 
order to suppress the resulting high frequency error introduced by the interpolations between stereographic patches. 

We accomplish this by modifying the evolution equation (|2.8|) as follows. In the code, (|2.8p is expressed in terms of 
a compactified radial coordinate x = r/(R + r), where R is an adjustable scale parameter and X + has finite coordinate 
value x — 1 . The evolution in retarded time u is carried out in terms of the variable $ = xJ, which is regular at X + . 
Then the evolution equation (|2.8[) takes the form 



dj (1 - * + $) =S, (2.9) 
where S represents the right hand side terms. We add angular dissipation to the it-evolution through the modification 

d u ( (1 - x)$, x + $ ] + e u h 3 d 2 Wd 2 ( (1 - x)$ >x + $ ) = S, (2.10) 



where h is the discretization size, e u > is an adjustable parameter independent of h and W is a positive weighting 
function with W = 1 inside the equator and W = at the patch boundary. This leads to 



d u (J(l - x)$ <x + $| 2 J + 2e u h 3 ${ ^(1 - x)$ iX + $J 9 2 W8 2 ^(1 - x)$ iX + $J } = 25R{ f (1 - x)* >x + $J S}. (2.11) 

Integration over the unit sphere with solid angle element dQ then gives 

d u j>\(l- x)$, x + <P\ 2 dil + 2e u h 3 j> W|8 2 ^(1 - x)$ lX + $^ \ 2 d£l = 23ft j> ^(1 - x )$ >x + $\ SdVi. (2.12) 

Thus the e^-term has the effect of damping high frequency noise as measured by the Li norm of (1 — x)^> x + $ over 
the sphere. 

Similarly, dissipation is introduced in the radial integration of (12. 9|) through the substitution 

d u (V - a:)$,x + du f (1 - x)$, a + + e x h 3 <5 2 WW<$> tU (2.13) 

with e x > . Angular dissipation is also introduced in the hypersurface equations through the substitutions 

(r 2 QX r -► {r 2 Q), r + e Q h 3 mWdBrQ (2.14) 

v r -> v r + e v h 3 mwmv. (2.15) 
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III. WAVEFORMS AT J+ 



For an analytic treatment of the Penrose compactification of an asymptotically flat space-time, it is simplest to 
introduce an inverse radial coordinate £ = 1/r, so that future null infinity X + is given by £ = [40l |. In the resulting 
x 1 * = (u, £, x A ) conformal Bondi coordinates, the physical space-time metric g^v has the conformal compactification 
9^1/ = £ 2 g^, where g pv is smooth at I + and, referring to (|2.ip . takes the form 

g^dx^dx" = - (e 2/3 V7 3 - h AB U A U B ) du 2 + 2e 2p dud£ - 2h AB U B dudx A + h AB dx A dx B . (3.1) 

The inverse conformal metric has the non vanishing components g ul — e~ 2 ® , g u = e~ 2/3 £ 3 V 7 g lA — e~ 2/3 U A and 

g AB = h AB 

The Bondi mass, news function and ty® (functions of u and x A ), which describe the total energy and radiation 
power, are constructed from the leading coefficients in an expansion of the metric in powers of £. The requirement of 
an asymptotically flat vacuum exterior imposes relations between these expansion coefficients. In the g pl , conformal 
frame, the vacuum gravitational equations are 

- £ 2 G^ = 2£{V^ V £ - g^V a VJ) + 3$ M „(V Q *)V^ (3.2) 

in terms of the Einstein tensor G^ v and covariant derivative V M associated with g^ v . Asymptotic flatness immediately 
implies that g u = (V a £)W a £ = 0(£) so that I + is null. From the trace of (|3.2p . we have 

(V a £)V a £ = -£Q + 0{£ 2 ), (3.3) 

where 

6 := V" V?t = e- 2f) ( d e (£ 3 V) + d A U A ) (3.4) 



is smooth at I + . In addition, (13. 2p implies the existence of a smooth trace-free field S M „ defined by 

£% v := V^VJ - jv,9. (3.5) 

For future reference we introduce an orthonormal null tetrad (n M , £**, m^) be such that n M = V M £ and £^d\x = di at 
J+. Note that (O, dill) and (JSU) imply 

m v mP{t vp + ^G up ) = 0. (3.6) 

The gravitational waveform depends on the value of £ MW on I + , which in turn depends on the leading terms up to 
0{£) in the expansion of g^. We thus expand 

h AB (u,e,x c ) = H AB (u,x c )+£c AB (u,x c )+0{£ 2 ). (3.7) 

Further conditions on the asymptotic expansion of the metric can be extracted from (|3.2j) . We have 

0(u,£,x c ) = H{u,x c ) + 0{£ 2 ) (3.8) 

(where the 0{£) term vanishes), 

U A = L A + 2£e 2H H AB D B H + 0(£ 2 ), (3.9) 

and 

£ 2 V = D A L A + £(e 2H K/2 + D A D A e 2H ) + 0(£ 2 ), (3.10) 

where 1Z and D A are the 2-dimensional curvature scalar and covariant derivative associated with h AB . These results 
combine with (|3.4p to give 

-2fl n tA i el t> i o„-2ff n riA„2H 



6 = 2e~ D A L +e[n + 2,e- zn D A D A e ZH + 0(t). (3.11) 
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In addition, the requirement that 

l{t AB - ^H AB H CD Z CD ) 

vanishes at T + implies via (|3.5[) that 

2H C(A D B) L C + d u H AB - H AB D C L C = 0(£). (3.12) 

The expansion coefficients H, H AB , c AB and L A (all functions of (u,x A )) completely determine the radiation 
field. One can further specialize the Bondi coordinates to be inertial at Z + , i.e. have Minkowski form, in which 
case H = L A = 0, H AB — q AB (the unit sphere metric) so that (|3.12p is trivially satisfied and the radiation field 
is determined by c AB . However, the characteristic extraction of the waveform is carried out in null coordinates 
determined by data on the inner worldtube so that this inertial simplification cannot be assumed. 



A. Calculation of the news 

The following calculation of the Bondi news streamlines the presentation in [f| and corrects errors. In order to carry 
out the calculation in the g^y computational frame, it is useful to refer to an inertial conformal Bondi frame [ifjj with 
metric g^y — fl 2 g^y = Lo 2 g^ v , where Q — u>£, which satisfies the gauge requirements that Q AB :— g AB \i+ = lo 2 H ab 
is intrinsically a unit sphere metric at 2 + and that (V a f2)V a f2 = 0(il 2 ). (See [4l| for a discussion of how the news 
in an arbitrary conformal frame is related to its expression in this inertial Bondi frame.) 

I + is a null hypersurface with the null vector h a — g a ^V/3f2|x+, or equivalently, h a — g a @V ' pl\x+ — ujn a , tangent 
to its generators. In order to complete a basis for tangent vectors to T + , let Q a be a complex field tangent to I + 
satisfying Q a n a — 0, g a pQ a Q l3 \i+ — and g a fiQ a QP\x+ = 2. In an inertial conformal Bondi frame, the news function 
can then be expressed as [j| 

N = lim ^-Q a Q (i V a Vf } n (3.13) 

evaluated in the limit of X + . (Our conventions are chosen so that the news reduces to Bondi's original expression in 
the axisymmetric case [l[). In terms of the g a p frame, with conformal factor i = 0,/uj, we then have 

N = \i m ±Q«QP(^LM - cuW a W p - + -L.g Q/3 (V^)V^ (3.14) 

= \Q a Q P (t a p-LoVjj [j - + -{dfg afi ){V^)V ll u J \. (3.15) 

A \ LO LO J 

(This corrects an error in equation (30) of 5j.) We determine w on I + in the g a p frame by solving the elliptic equation 
governing the conformal transformation of the curvature scalar (|2.7j) of the geometry intrinsic to &u — constant cross- 
section to a unit sphere geometry, 

K = 2(uj 2 + H AB D A D B logu). (3.16) 
The condition that (V a f2)V Q f2 = 0(f2 2 ) determines the time dependence of u, 

2h a d a logw = -e- 2H D A L A , (3.17) 

which is used to evolve u> given a solution of (|3.16p as initial condition. 

In order to obtain an explicit expression for the news (|3 . 1 5|) in the g a ^ frame we need to fix the choice of . The 
freedom Q@ — > Q? + \n^ leaves (|3. 15|) invariant but it is important for physical interpretation to choose the spin 
rotation freedom QP — > e~ la Q? to satisfy n a \7 a Q 13 = O(O), so that the polarization frame is parallel propagated 
along the generators of T + . This fixes the polarization modes determined by the real and imaginary parts of the news 
to correspond to those of inertial observers at X + . 

We accomplish this by introducing the dyad decomposition H AB — (F A F B + F A F B )/2 where 




F A = Q \j - Q J\ oTttVtt- (3-18) 



7 



We set Q? = e~ lS uj~ 1 F 13 + Xh 13 , where F a := (0, 0, F A ). The requirement of an inertial polarization frame, h a V a Q 13 = 
0(Q), then determines the time dependence of the phase S. We obtain, after using (|3.17[) to eliminate the time 
derivative of oj, 

2i(d u + L A d A )5 = D A L A + H AC F c ((d u + L B d B )F A - F B d B L A ). (3.19) 
We can now express the inertial news (|3 . 1 5|) in the g a p frame as 

N = l e - 2i6 uj- 2 F a F' 3 (t af3 - uV a V p - + -(d e g a0 )(W^)W^uA . (3.20) 

Z y UJ UJ J 

with F a — (0, 0, F A ). An explicit calculation leads to 

N = - A e-^uj- 2 e- 2H F A F B {(d u + £ L )c AB - \c AB D c L c + 2ujB < A {u;- 2 D B (ue 2H )}}, (3.21) 

where £l denotes the Lie derivative with respect to L A . This corrects a minus sign error in (38) of 5], where 
spin-weighted expressions for the terms in (|3.21[) are given. 

In inertial Bondi coordinates, the expression for the news function reduces to the simple form 

N =\Q A Q B d u c AB . (3.22) 

However, the general form (|3.21|) must be used in the computational coordinates, which is challenging for maintaining 
accuracy because of the appearance of second angular derivatives of uj. 



B. Calculation of the Weyl tensor 

Asymptotic flatness implies that the Weyl tensor vanishes at X + , i.e. Cp, vpa = 0(1) in the g^ v conformal Bondi 
frame (13. 1|) . This is the conformal space version of the peeling property of asymptotically flat spacetimes [3J|. In terms 
of the orthonormal null tetrad (n M , £^ , m M ) , with n M = and i^dfi = di at I + , the radiation is described by the 
limit 



* := -\ lim \r^m v ffnfC^ p<T , (3.23) 
2 e^o £ * 1 



which corresponds in Newman-Penrose notation to —(1/2)^. The limit is independent of how the tetrad is extended 
off T + but to simplify the calculation we make the following choices adapted to our conformal Bondi coordinates. We 
set & = e 20 V p u, = + 0(1), l p V p rh^ = and I^V p n p = 0, which implies 

% = Vpl- ^£ P ® + 0(e 2 ). (3.24) 

Our main calculational result is: 

* = \wm u mp(v»± up -V v ±^\ x+ , (3.25) 
and that 1|3.25[) is independent of the freedom 

rn v — > rh v + \n u . (3.26) 

The result f|3 . 25[) follows from the following sequence of calculations beginning with (13 . 23[) (where evaluation at X + 
is assumed): 

-2* = 1 n"m u h p m a C flvp(T (3.27) 



1 

e 

i 



h^m"h p m' T R^ p<7 (3.28) 



Wm v rh, p (V ll V v fi p -V v V M fi p ) (3.29) 
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1 / - - - » £<3 - - £& 

-h»rh v m? V - V v (ffi w ) - V M V,(— £ p ) + V„V M (— * p ) ) (3.30) 



w 

- -n'WWM E^-S^-V^V^j ^) + V„V„(— *„) ) (3.31) 

- ^"m"m' ( - E w ?t - V.if^) + V„V„(*/ P ) | (3.82) 
_ ^Wft^-n"^^ (3.33) 



m"W ^S„ p - hWti^&jp - (3.34) 



jm^fs^ + ic^) (3.35) 



= -ft"mW^ r V»E w j. (3.36) 

Here (I3.28P follows because all trace terms vanish; p. 291) follows from the commutator of covariant derivatives; (|3.30|) 
follows from (|3 . 5p : (|3.31|) follows from differentiation; (|3.32p follows from ()3.24[) : (|3.33p follows from taking leading 
terms and using the covariant commutator; (|3.34[) follows from the vanishing of the Weyl tensor at T + ; (|3 . 35[) follows 
algebraically; and (|3.36[) follows from (|3.6p . 

Invariance of ^ under the freedom (|3.26p follows from noting that 

n p 'm v h p h a C lil/pa = (3.37) 
and then following the steps analogous to (|3.27p - (|3.36p to show 



n»rh v n p ^V M £„„ - V,Sp P J \ I+ = 0. (3.38) 

Finally, the Weyl tensor must be scaled intrinsic to the g pv conformal frame in order to describe the radiation 
observed by inertial observers at T + . The conformal transformation g pv — ui 2 g pv gives for the inertial radiation field 

* : = -\ lim h^Q v nPQ a C pvpa 

= -^u- 3 e- 2iS Kmjn»F v nPF'>C fiVp(T , (3.39) 

where Q? = e- lS uj- 1 F 13 + Xn 13 is the same inertial polarization dyad used in describing the news p. 211) . From (|3.25p . 
we then have 



* = ^ui~ 3 e~ 2iS n^F u F p (vjl vp - V,E^ \ x + . (3.40) 
We next need to express in terms of the computational variables. The straightforward way is to expand (|3.40p as 
* = ^- 3 e- 2tS n^F A F B (d p t AB - d A t pB - T a pB t Aa + \ I+ . (3.41) 
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and calculate the individual components of S M „ in terms of those variables. This involves lengthy algebra, which is 
simplified by the following intermediate results which hold at X + : 

t u = -2d £ 2 /? (3.42) 
Em - \e- 2H d t {h AB d t U B ) (3.43) 

t tu = \e 2H Tl + ^D A D A e 2H - L A ± eA (3.44) 

(V^XV^XV = \er 2H {d u + L A d A )(e- 2H D A L A ) (3.45) 

(V^)S M = l -d A {e- 2H D B L B ) (3.46) 

±ab = \e- 2H (d u + C L )c AB + e- 2H D A D B e 2H - jH AB (K + Zer 2H D c D c e 2H ). (3.47) 

We use a Maple script to convert these expressions in terms of 3 operators acting on the spin- weighted computational 
fields and construct the final Fortran expression for ^ . 
In inertial Bondi coordinates, (|3.4ip reduces to 



* = \Q A Q B dlc AB = d 2 u d t J\ + . (3. 



This is related to the expression for the news function in inertial Bondi coordinates by 

* = d u N. (3.49) 

However, as in the case of the news, the full expression (I3.41[) for must be used in the code. This introduces 
additional challenges to numerical accuracy due to the large number of terms and the appearance of third angular 
derivatives. 

C. Linearized expressions 

The general nonlinear representation of \& in (|3.4ip in terms of the computational variables is quite long but reduces 
to a simpler form in the linearized approximation, i.e. to first order in perturbations off a Minkowski background. In 
terms of the spin-weighted fields J = h AB q A Q B /2 and L = L A q A , we find 

* = ~t tu BL + d u tj - ^dt u - X -d u J (t lu + t K ) (3.50) 

(evaluated at X + ) , where the only nonvanishing zeroth order parts of £^„ are 

= \q A Q B ^AB = -\ %u = \ (3.51) 

and the required first order components are 

tj = ^q A q B ± AB = S 2 H - i J + X -3 u d t J (3.52) 

t u = q A t uA = ig 2 Z + -mL + \l. (3.53) 
4 4 2 
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Then (|330| reduces to 

* = \dld t j - ^d u j - igx - ig 2 (9Z + Bl) + d u &H. (3.54) 

In the same approximation, the news function is given by 

N = p A q B (d u c AB + 2D A D B {u + 2H) \ = ^d u d e J + hi\uj + 2H). (3.55) 
Using the asymptotic relations 

d u J = -dL (3.56) 

d u uj = -i(9L + 8L), (3.57) 

which arise from the linearized versions of (|3.12p and (|3.17p . it is easy to see that (I3.49|) . i.e. W = d u N, still holds in the 
linearized approximation. (In the nonlinear case, the derivative along the generators of X + is n^d^ — e~ 2H (d u + L A d 'a) 
and (|3.49p must be modified accordingly.) 

The linearized expressions (13.54[) and (|3 . 55[) provide a starting point to compare the advantages between computing 
the radiation via the Weyl component '5 or the news function TV. The troublesome terms involve L, H and lo, which 
all vanish in inertial Bondi coordinates. One main difference is that ^ contains third order angular derivatives, e.g. 
<3 3 L, as opposed to second angular derivatives for N . This means that the smoothness of the numerical error is 
more crucial in the \& approach. Balancing this, another main difference is that N contains the d 2 ui term, which is a 
potential source of numerical error since u must be propagated across patch boundaries via (|3 . 1 7|) . 



D. Summary of the gravitational radiation calculation 



The characteristic Einstein equations are evolved in a domain between an inner radial boundary at the interior 
worldtube, and an outer boundary at future null infinity. Initial data for J(u,r, x A ) is required at u = 0. This data 
is constraint-free so that, in the absence of an exact solution or other prescription of the data, we can simply set 
J(0,r, x A ) — 0. Alternatively, in order to reduce spurious initial radiation, we can set the Newman-Penrose Weyl 
tensor component ^(O, r, x A ) = 0, which determines J(0,r,x A ) = when continuity conditions are imposed at the 
inner worldtube. The metric data from a Cauchy evolution is interpolated onto the inner worldtube to extract the 
boundary data for the characteristic evolution. This extraction process involves carrying out the complicated Jacobian 
transformation between the Cartesian coordinates used in the Cauchy evolution and the spherical null coordinates 
used in the characteristic evolution. The full details are given in [l(j. The result is boundary data for J, f3,U,Q,V 
on the worldtube, which supply the integration constants for a radial numerical integration of (|2.3p , (|2.4[) . (|2.5[) and 
(|2.6p . in that order. Given the initial data J(0,r, x A ), this leads to complete knowledge of the metric on the initial 
null cone. Then (|2.8p gives an expression for J ur , which is used to determine J on the "next" null cone, so that the 
process can be repeated to yield the complete metric throughout the domain, which extends to T + . 

Before the gravitational radiation is calculated from the metric in the neighborhood of X + , it is necessary to compute 
the auxiliary variables oj(u, x A ) and 6{u, x A ) which determine the inertial polarization dyad in which to measure the 
news function N or Weyl component 'P. Given a solution of (|3 . 1 6|) for the initial value of ui(0,x A ), its evolution is 
computed by integrating (|3.17p . (If J = initially, then uj = 1 is the solution to (|3.16p . Otherwise, oj is initiated by 
solving a 2-dimensional elliptic equation.) Similarly, fixing the initial polarization basis by 6(0, x A ) — 0, its evolution 
is computed by integrating (|3.19p . Then the news N is given by (|3.2ip (or in spin- weighted form by the formulas in 
Appendix B of Ref . Q ) and \& is given by (|3.4ip . 

The above procedure computes TV or ^ as functions of the code coordinates (u, x ), rather than inertial coordinates. 
In the linearized case, which is used for the tests in Sec. IVI1 the change to inertial coordinates is a second-order effect 
that can be neglected. However, that is not the case in general and the full procedure is described in Sec. IV B of 
Ref. @. 



IV. PATCHING THE SPHERE 



The nonsingular description of smooth tensor fields on the sphere requires more than a single coordinate patch. 
Here we consider the stereographic treatment which uses 2 coordinate patches, and the cubed-sphere treatment, which 
uses 6 patches. In both cases the metric qAB of the unit sphere is expressed in terms of a complex dyad q A (satisfying 
q A QA = 0, q A q~A — 2, q A — q AB qs, with q AB qsc = S A and qAB = \ [qAqs + qAqs))- The dyads for each patch are 
related by spin transformations at points common to more than one patch. 
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A. Circular stereographic patches 



In stereographic coordinates, the sphere is covered with two patches, one for each hemisphere. The North hemisphere 
is covered by the complex stereographic coordinate £jy = Vn + Wn, which is related to standard [6, <f) angular 
coordinates by £jy = tan(6>/2)e*^ and which is regular on the entire sphere except for the South pole. The South 
hemisphere is covered by the complex stereographic coordinate £s = 1/£n — Vs + ips = cot(0/2)e _t ^, which is 
singular at the North pole. Every point on the sphere is covered by at least one of the patches, and there is a region 
around the equator where points are covered by both patches. In this overlap region between the two patches, a scalar 
function F with value Fn(^n) on the North patch has the value Fs(£s — l/6v) on the South patch. For a function 
F of spin-weight a, F s (^ s = l/6v) - F N (Z N )(-l) s e- 2is * . 

In the x A — (77, p) coordinates, the unit sphere metric in each patch is given by 

q AB dx A dx B = -^(d V 2 + dp 2 ), (4.1) 



where 



The equator corresponds to the circle 



We fix the dyad by the explicit choice 



P = l+r] z +p z . (4.2) 



VV+7=1- (4-3) 



9 A = f(M), » = V=1. (4.4) 

In the composite mesh method, all boundary points of one patch are interior points of another patch. The over- 
lapping of the patches is key to the stability of this method. The two stereographic coordinate patches must both 
extend beyond the equator. In the scheme originally used to implement the computational S-formalism j30l | in the 
characteristic code, the North and South patches were represented by square (77, p) grids. In the scheme implemented 
for meteorological studies [28[, circular masks are applied so that the computational grids extend only a few zones 
beyond the equator. Here we adopt this circular grid boundary but place it a fixed geometrical distance past the 
equator, i.e. the grid boundary for the North patch is a circle lying in the South patch. The finite buffer zone 
between the equator and the grid boundary allows for angular dissipation, as developed in Sec. Ill Al to damp the high 
frequency intergrid interpolation error before it crosses the equator. This protects measurements of the news function 
(or ty) in the North patch, which involve two (or three) angular derivatives, from substantial contamination by the 
interpolation error at the patch boundary. 

We discretize the stereographic coordinates according to 

m = -l + (i-0-l)A, Pj = -l + (j-0-l)A (4.5) 

where, following the notation in [28| . O is the number of points (overlapping points) by which each grid extends 
beyond the equator and the indices range over l<i,j<M + l + 20, with M 2 being the number of grid points inside 
the equator. The grid spacing A depends on M according to 

A=A (46) 

The square grid determined by (|4.5|) ranges over 

( Vl , P] )e(-l-OA,l + OA) (4.7) 

in each patch. 

In the original square patch method (30j , the evolution algorithm is applied to the entire set of points in the square 
(77, p) grid, with the field values at the resulting ghost points supplied by interpolation from the other patch. In the 
circular patch method 28], the evolution algorithm is only applied to points inside a circle r = y/r/ 2 + p 2 , where 
r > 1 so that the boundary lies a small distance past the equator. In convergence tests, the number of overlap points 
determined by O is adjusted so that r is at a fixed position for all grids, i.e. O scales as 1/A. The grid points outside 
this circle are either ghost points or inactive. The circular patch is clearly more economical than using a square patch 
and avoids the error introduced by the large stereographic grid stretching near the corners of the square. 
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When the finite-difference stencil is used near the boundary of the active grid points, field values required at the 
ghost points outside a circular patch are interpolated from values at interior points of the opposite patch. The 
algorithm for determining the value of a scalar function F/v at a ghost point in the North patch starts with the 
determination of the ghost point's coordinates in the overlapping South patch, followed by the interpolation of the 
value of the function Fs at the ghost point, i.e. the Fn ghost point values are obtained by interpolation via the Fs 
active grid values. 

Let Re be the width of the finite-difference stencil divided by 2A. In the circular patch method, we define the 
active finite difference grid, i.e. the grid points to which the evolution algorithm is applied, by 

yj'fi + P] < 1 + (O - Re) A, (4.8) 

where O > Re- Stability of the composite mesh method requires that the interpolation stencil for the ghost points 
for one patch lies below the equator in the other patch. Those requirements give a minimum value of O but a larger 
value may be necessary to establish an effective buffer zone for the dissipation to attenuate the interpolation error 
before it enters the opposite patch. The optimum value of O in order to avoid instability or inaccuracy, needs to be 
established by experiment. (Too large a value would lead to inaccuracy due to the stretching of the stereographic 
grid.) 



B. The cubed sphere 

In the cubed-sphere approach, developed for meteorological studies in (29| and later for numerical relativity in [3l| , 
6 coordinate patches on the sphere are obtained by projecting the 6 faces of a circumscribed cube. The method has 
recently been applied to characteristic evolution in [33 ] and independently in [33| . Here we follow the notation of [32| , 
except we denote the angular coordinates by (0i,02) (rather than by (/?, er)). In addition, in order to ensure that 
the coordinates and dyads on each patch are consistently right-handed, with the vector cross-product vector pointing 
out of the sphere, we introduce some sign changes in the conventions used in [32| for the coordinate transformations 
between the patches and in the definition of the dyad q A . These conventions simplify the interpatch transformation 
of spin-weighted quantities. 

Given Cartesian coordinates (x, y, z), we define angular coordinates x A = (0i, 02) on the 2-sphere x 2 + y 2 + z 2 — 1 
by means of the six patches (x±, y±,z±), where 



x± : 0i = arctan , 02 = arctan ^— J 

y± : 0i = arctan [ ±— ) , 0a = arctan ( — 



z± : 4>\ = arctan > 4>2 = arctan ^— ^ . (4.9) 

In each patch, the range of the coordinates is — 7r/4 < 0i, 02 < tt/4 and the metric is 



ds 2 = (l - sin 2 0i sin 2 2 ) 2 ^cos 2 02 #? + cos 2 0i dxfe - ^ sin( 2o, i sin( 2'> 2 ) ,{,,., ) . ( 4. 1.0) 

As a simple dyad representing (14. 1 0[) . we choose 

id s ) cos 02 (6> c -H6> s )cos0i \ A ( na n 6 S 
HA ~\ 49 2 e 2 ' 49 2 9 2 )> r 

where 



Z202 > Imffl h 1 ~ 2^ s ^-p,20A^-p , (4.11) 
V w£vj wivi J V cos 02 cos 0i ' 



1 — sin 0i sin 02 / 1 + sin 0i sin 02 



(4.12) 



The operator 3 acting on a field S with spin-weight s is (5S = q A dAS + sTS where, with the present conventions, 

cos 2 0i cos 2 2 (sin 0i + sin 2 ) + (cos 2 0i — cos 2 02 ) (sin 2 — sin 0i ) 



40 c COS 02 cos 01 

. cos 2 0i cos 2 02 (sin 0i — sin 2 ) + (cos 2 02 — cos 2 0i ) (sin 0i + sin 2 ) 

A9 S COS 02 COS 01 



(4.13) 
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We introduce ghost zones in the usual manner along the boundaries of each patch, and couple the patches together 
by interpolating the field variables from neighboring patches to each ghost point. With the definition (|4.9[) . the 
angular coordinate 4>\ or <j>2 perpendicular to an interpatch boundary is always common to both adjacent patches. 
This greatly simplifies interpatch interpolation, since it only needs to be done in 1 dimension, parallel to the boundary. 



V. COMPARISON BETWEEN STEREOGRAPHIC AND CUBED-SPHERE METHODS 

We carry out a test of wave propagation on the sphere to compare the accuracy of using circular stereographic 
patches with the cubed-sphere methods, with emphasis on the accuracy of the angular derivatives required in waveform 
extraction by the news and ^4 approaches. The test allows direct comparison between the stereographic and cubed- 
sphere treatments without introducing the complications of characteristic evolution and the explicit calculations of 
the news or ^ . 

The test is based upon solutions to the 2D wave equation 

-<9 2 $ + 99$ = 0, (5.1) 



where $ = cos{ujt)Yi m , ui — \J i(t + 1) and Yi m are spherical harmonics. 

For the case I = m = 2, we compare test results for the stereographic grid with circular patches and the cubed- 
sphere grid. For the stereographic grid, the simulations are run with M 2 grid points in each patch, for M — 100 and 
M = 120. The corresponding cubed sphere runs keep the number of grid-cells covering the sphere the same as for the 
stereographic case, not counting those cells that overlap with another patch. For M 2 stereographic grid points there 
are ~ irM 2 /4 grid-cells inside the equator on each hemisphere. In the cubed sphere grid, with iV 2 points per patch, 
the entire sphere is covered by 6 x TV 2 points. Equating the number of cells for the entire sphere gives N 2 ss (7r/12)M 2 . 
The above values of M then correspond to N = 51, 61. The tests are run until t = 120. 

Angular dissipation is necessary for the stability of the stereographic runs. For grid size A, it was added in the 
finite-difference form 

<9 2 $ -> d 2 <f> + eA 3 X> 4 9 t $, (5.2) 



where L> 4 $ = ^-^-(L> +r) P_ ?) + 2? +p L>_ p )^ $, where V + (or T> ) indicates the forward (or backward) difference 

operator in the indicated direction. Experimentation with tuning the dissipation revealed that a small value e = 0.01 
is sufficient to suppress high frequency error. The finite difference stencil (taking dissipation into account) has width 
Re = 2. Using a Ath order Lagrange interpolator and by tuning the number O of overlapping points, we obtained 
good results with O = 5. Angular dissipation was not used in the cubed sphere runs. 
We use the L m norm to measure the error 

£{Q) = I \&numeric ^analytic \ | oc ■ (^-3) 

We measure the convergence rate for £ ($) at a given time t, for two grid sizes Ai and A2, by 

_ log 2 (£(<fr) Al /g(<fr) Aa ) 

log^AVAx) ' (5 - 4) 

Convergence rates for other quantities are measured analogously. For a given grid, we measure the error for the circular 
patches in the North hemisphere; while for the cubed-sphere method we measure the error on the (+£, +y, +z) patch, 
excluding ghost points at the edges of the patch. The finite difference approximations for the codes are designed to 
be second order accurate. 

Excellent second-order convergence of £ ($), based upon the M — 100 and M = 120 grids, is evident for both 
methods from the results listed in Table Q] The time dependence of the error plots in Fig. [IJ shows that the cubed- 
sphere error is ~ ^ the stereographic error. 

A more important test to assess the error relevant to gravitational wave extraction is to measure the error in 
S 2 $, since second angular derivatives enter in the computation of the Bondi news. The convergence rates, measured 
with the Lqo norm, are shown in Table [XTJ The circular patch and cubed-sphere results indicate clean second order 
convergence up to the final run time at t = 120. The plots of the error versus time in Fig. [2] show that the error for 
the cubed-sphere is about |rd the stereographic error. 

Similarly, it is important for the purpose of gravitational wave extraction using the Weyl tensor to measure the 
error in 99 2 <i>, since third angular derivatives enter into the computation of \&. The corresponding convergence rates 
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ALGORITHM 


t=1.2 


t=12 


t = 102 


t = 120 


circular patch 


2.002 


1.988 


1.994 


1.999 


cubed-sphere 


1.994 


1.970 


1.982 


1.985 



TABLE I: Convergence rates for £($), obtained with the norm using the M — 100 and M = 120 grids. 




ALGORITHM 


t=1.2 


t=12 


t = 102 


t = 120 


circular patch 


2.022 


1.945 


1.992 


2.006 


cubed-sphere 


1.954 


2.019 


1.997 


1.971 



TABLE II: Convergence rates for the Loo error £ (9 <E>) for various times t, obtained using the two highest resolution runs. 




FIG. 2: The Loo error £(3 2< I?) vs t for the highest resolution runs is compared using the circular patch method and the 
cubed-sphere method. 
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are shown in Table IIIII Now the cubed-sphere method shows poor convergence of the L m error at early times. The 
underlying error is generated at the corners where three patches meet, as indicated by the improved convergence rate 
measured using the Li norm. Apparently some built-in dissipation of the evolution algorithm smooths this patch- 
boundary error and second order convergence is evident by by t = 102. To a much smaller extent, the L x error for 
the circular patch method also shows some deviation from second order convergence at early times, but clean second 
order convergence is evident by t — 12. 

The magnitude of the L x error in §3 2 $ is plotted vs time in Fig. [3] Until about t = 60, the cubed-sphere method 
has the largest error. But at the end of the run at t = 120 the cubed-sphere error is about | the stereographic error. 

Surface plots of the error at the final run time are shown in Fig. [4] The circular patch and cubed-sphere errors are 
both smooth, as would be expected of the second order truncation error arising from the finite differencing. For the 
circular patch, this shows that dissipation in the buffer zone surrounding the equator effectively guards against the 
high frequency error introduced at the patch boundary. Our results for the stereographic method justify its use in 
the comparison of the news and Weyl tensor extraction strategies in Sec. IVII 



ALGORITHM 


t=1.2 


t=12 


t = 102 


t = 120 


circular patch, Loo norm 


2.278 


2.032 


1.988 


2.009 


cubed-sphere, Loo norm 


1.108 


0.882 


2.009 


1.959 


cubed-sphere, L2 norm 


1.883 


1.983 


1.981 


1.959 



TABLE III: Convergence rates for the error £ (33 2 $). For the cubed-sphere method the dominant error arises at the patch 
corners, which is revealed by the comparison of the Loo and L2 errors. This can be understood in terms of the the inter- 
patch interpolation stencil which is partially off-centered near the corners, where the error is greatest. The inherent numerical 
dissipation of the evolution algorithm keeps this localized, non-smooth noise from growing, while smoother error from other 
regions grows linearly in time (see Fig. [3} . The net effect is that at late times both the Loo and the L2 norms of the error in 
the third derivative show second order convergence, while early in the evolution the Loo shows only first order convergence. 




40 80 



t 

FIG. 3: The error £(33 2( I>) vs t is compared for the highest resolution runs using the circular patch method and the cubed-sphere 
method. 




FIG. 4: 2D snapshots of the error in 93 2 $ at t = 120 for the on the North hemisphere for the circular patch method (left plot), 
and for a cubed-sphere patch (right plot). For the sake of plot clarity these 2D snapshots use only every third data-point along 
each axis. The third angular derivatives are smooth for both methods 
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VI. COMPARISON OF NEWS AND WEYL TENSOR EXTRACTION 



Here we compare the accuracy of waveform extraction by computing the news function N or the Weyl tensor 
component Sfr in a linearized gravitational wave test problem. The computations are carried out by the procedure 
described in Sec. MI Dl In accord with (|3.49[) . the \& computation yields an alternative numerical value for the news 

N*=N\ U=0 + f Vdu, (6.1) 
Jo 

where N — in the analytic problem. We compare the two extraction methods in terms of the errors in N and 
obtained using the stereographic method. 

We base the test on a class of solutions in Bondi-Sachs form to the linearized vacuum Einstein equation on a 
Minkowski background given in Sec. 4.3 of [42j. The solution allows us to make convergence checks of the Bondi- 
Sachs metric quantities as well as the news function. The solutions are expressed in terms of spin-weighted spherical 
harmonics s Ye m [H, [43[, modified to avoid mixing of the to and — m components when extracting the real part 
according to [lij 

1 i 
s R tm = -j= [ s Y lm + (-l) m s Y^_ m ] for to > 0, s R lm = -j= [{-l) m s Y im - s Y^ m ] for to < 0, s R m = s Y m . (6.2) 

Ref. [H| gives explicit expressions for the s Re m in stereographic coordinates. 
Following [32], these linearized solutions have Bondi-Sachs variables 



J = y/(i - 1)1(1 + !)(£ + 2) 2 R em ^(J e (r)e l " u ), U = iR im ^{Ui{r)e ivu ), 

= R lr Jt{0 t e ivu ), W c = R lm sR(W ci (ry m ), (6.3) 

where W c determines the perturbation in V. Here Ji(r), Ui(r), 0e, W c g(r) are in general complex, and taking the real 
part leads to cos(z^u) and sin(^it) terms. The quantities and W c are real, while J and U are complex. We require a 
solution that is well-behaved at future null infinity, and is well-defined for r > ro > 0, where r$ is the inner boundary. 
We find in the case £ = 2 



Mr) 
U 2 (r) 



02 = 0o 

24/3 + 3iuCi - iv z C 2 C x C 2 



36 4r 12r 3 

-2tev0 a + 3v 2 G\ - v 4 C 2 20 o C x ivC 2 C 2 



36 r 2r 2 3r 3 4r 4 

24iu0 o - 3^ 2 Ci + iy*C 2 MuC x - 60 o - iv 3 C 2 v 2 C 2 ivC 2 C 2 



6 3r r 2 r 3 2r 4 ' 

with the (complex) constants 0o, G\ and C 2 freely specifiable. In the case I = 3 

02 = 00 



Mr) 
Mr) 



2 



6O0o + 3wCi + v^C 2 Cx ivC 2 Q 
180 + Wr ~ 6r 3 ~ Ar 1 

-6Ow0o + 3is 2 d - iv 5 C 2 20q C x 2v 2 C 2 bwC 2 C 2 



180 r 2r 2 3r 3 4r 4 r 5 

6Oiv0 o - iv 2 C x + iv 5 C 2 ivCx - 20 o + v A C 2 i2v i C 2 Mv 2 C 2 5vC 2 3C 2 



WMr) = ^ + ^ ~ 2 jsr- + — + — ■ (6.5) 

The news TV for the linearized wave is given by 



£(£ + 1) 



V - n ( ( Hm [ - ~ v " - ^ J t - Y r2 Mj + e iuu 0ej y/(e -!)£(£+ !)(£ + 2) 2 R em , (6.6) 



corresponding For the cases I —2 and 3, this gives 



Y ■)," ( %J ^=-e w A 2 R 2m - £ = 3: N = U^ e lvu ) , /?,,„ . ( (i.7 ) 



For the linearized case ^ = 7V„, which gives 



I = 2 : * = . )f ( e i™\ 2i?2m . t = 3 . ^ = R ^ ^ 2 e wu j ., /?:„„ . ( I..S ) 
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A. Test specifications 

Tests were run with the solution parameters v = 1 and m = for the cases I = 2 and I = 3, with 

Ci = 3-l(T 6 , C 2 = 1(T 6 , /3 = MCr 6 (* = 2) (6.9) 
Ci=3-10" 6 , C 2 =i-10~ 6 , /3 = i-10~ 6 (i = 3). (6.10) 

The inner worldtube boundary was placed at ro = 2 corresponding to a compactified radial coordinate x$ = tq /{R + 
ro) ~ .1888, where we have set the scale parameter R = 9. 

For the convergence measurements, the (r),p,x) grid consisted of M 3 points, with M = 100 and M = 120. The 
boundary of the circular patches were fixed at yjrf + p 2 = 1.4. The runs were stopped at t = 100. The L m and L 2 
error norms were computed on the North hemisphere, using the values from the North patch. 

Angular dissipation was applied only to the circular patch runs, with the dissipation coefficients e x — 0.009, 
e u = 0.0009, €q = €w — 0.00001. The weighting function W for application of the dissipation was taken to be a unit 
step function which vanishes for \/r] 2 + p 2 > 1.3. 

We present output data for the real parts of J, N and Ny, For the m = case, these quantities correspond to a 
pure © polarization mode. For comparison purposes, we include results for the circular patch without dissipation and 
the original square patch treatment. 

B. Test results for J 

We first present test results for J, which is a typical metric quantity entering into the waveform calculation. Figure 
[5] show the Loo norm over the North hemisphere of the error £(J) vs the compactified radial coordinate x at the end 
of the run at t = 100 for the I = 2 wave. The figure compares runs made with the circular patch method (dissipation 
applied) with runs without dissipation and runs with the original square patch method. The plots show that angular 
dissipation reduces the error. This will become more evident in the later test results for the news in which higher 
angular derivatives are involved. An important feature of the plots is that in all cases the error increases monotonically 
and takes it maximum value at T + (x = 1), as would be expected of the radial marching algorithm. This allows us to 
focus our error analysis on output at X + . 

Table HVl gives the convergence rate of the error in J measured at T + at various times during the £ = 2 run for the 
three methods shown in Fig. [5] Clean second order convergence, measured either with an i 2 or Loo norm, is indicated 
in all cases. The corresponding convergence rates for the I = 3 runs are given in Table fVl The I = 2 runs are more 
discriminating because |J| has a sin 2 9 dependence which peaks at the equator close to the interpatch interpolation, 
as opposed to the sin 2 9 cos 9 dependence of the £ = 3 case which vanishes at the equator. In the following we restrict 
our discussion to the 1 = 2 case. 

The time dependence of the L 2 and Loo errors in J at X + for the circular patch run (with dissipation) is plotted 
in Fig. [S] The plots are based upon output at integer values of t, which samples the error at various phases during 
the underlying period T = 2n. The errors for the two grids used in the convergence measurements are rescaled to the 
values for the finest M = 120 grid, with the overlap again confirming clean convergence. The magnitude of the error 
is approximately 0.1% the value of J. The L 2 error is smaller than the Loo because the error is sharply peaked near 
the equator. This error pattern in the North hemisphere is exhibited in the snapshot of Fig. [7] at t — 100. The profile 
is quite smooth - some of the apparent jaggedness near the edge is an artificial effect of the irregular pattern of grid 
points at the edge of the equator. The sharp spikes in the corresponding error snapshot for the circular run without 
dissipation shown in Figs. [5] illustrate the essential role of angular dissipation in guarding the Northern hemisphere 
from the interpolation error at the circular patch boundary. Such spikes are not apparent in the corresponding error 
snapshot for the square patch shown in Fig. [5] The more regular square patch boundary does not require angular 
dissipation, although the resulting error is larger than for the circular patch with dissipation. 
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FIG. 5: The Loo error £(J) plotted vs x at t = 100 for runs with the circular patch method (with and without dissipation) and 
with the square patch method. 
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FIG. 6: Plots of the error £ (J) at X + vs t, for the circular patch with dissipation measured with the L2 norm (left plot) and 
the Loo norm (right plot) . The error for the M = 100 grid is rescaled and overlaid on the error for the M = 120 grid to exhibit 
the second order convergence. The smaller L2 error indicates that the maximum error arises near the equator. 



Variable 


circular patch 


circular, no dissipation 


square patch 


£i 2 (J)t=i 


2.01 


2.00 


2.01 


£l 2 (J)t=io 


2.01 


1.98 


2.00 


£l 2 (J)t=90 


2.00 


2.02 


2.02 


£l 2 (J)t=10Q 


1.92 


2.03 


2.00 




2.01 


2.01 


2.01 


£Loc {J)t=10 


1.95 


2.00 


1.99 


£L^{J)t=90 


2.07 


1.96 


2.00 


£Loc {J)t=10Q 


1.92 


2.01 


1.99 



TABLE IV: Convergence rates of the error £(J) at X + for the £ = 2 run, measured at times t = 1, t = 10, t = 90, and t = 100. 
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Variable 


circular patch 


circular, no dissipation 


square patch 


£L 2 (J)t=l 


2.02 


2.01 


2.01 


£l 2 {J)t=io 


2.00 


2.00 


2.00 


Sl 2 (J)t=90 


2.03 


2.02 


2.02 


£h 2 (J)t=WQ 


2.05 


2.00 


2.01 


£i M {J)t=l 


2.02 


2.02 


2.02 


£ioo {J)t = 10 


1.99 


1.99 


2.00 


£-L<x (J)t=90 


2.02 


2.02 


2.04 


£koo (J)t=10O 


2.00 


2.00 


1.99 



TABLE V: Convergence rates of the error £{J) at 1 + for the I — 3 run, measured at times t = 1, t = 10, t = 90, and t = 100. 




FIG. 7: Surface plot of the error in J at t = 100 for the circular patch run (with dissipation) . 



FIG. 8: Surface plot of the error in J at t = 100 for the circular patch run without dissipation. 




FIG 



9: Surface plot of the error in J at t = 100 for the square patch run. 



22 



C. Test results for the news 

We now compare test results for the news function in terms of a direct calculation of N via (|3 . 2 1 [) and a calculation 
of via (|6.ip using the Weyl component W given in (13.41[) . We restrict the discussion to the t = 2 runs which are 
more challenging than I = 3 with respect to problems near the equator. Tables ED and EED give the convergence rates 
of the L2 and errors in N and N-q, measured at various times for runs with the circular patch (with dissipation) , the 
circular patch without dissipation and the original square patch methods. At the final run time t = 100, measurements 
for all cases show clean second order convergence, although there is a small departure in the JV* rates at early times. 

The plots of the L2 error vs time for the circular patch runs in Fig. [TU] show little difference in the time behavior 
between N and although the error in is slightly smaller. The errors measured at the end of the runs on 
the M = 120 grid are given in Table rvTTTl for the circular patch, the circular patch without dissipation and the square 
patch. The best results are obtained for the circular patch, which shows an w 30% improvement over the original 
square patch treatment. The results also show the essential improvement due to the use of angular dissipation. For 
the circular patch, the error in was « 24% smaller than the error in N at the end of the run. But it is also clear 
from the plots of the L 2 error Fig. [Til] that this ratio depends when and where this ratio is taken. At the equator 
where the news takes its maximum value, its modulus for this test is | N ana i ytic | w 8 x 10~ 8 . At the end of run, the 
corresponding fractional errors in and N are « 4% for averaged values and rs 9% for the maximum errors at the 
equator. 

Surface plots of the error in N and A* at the end of the run are given in Figs. fTTI- 1141 for the circular and square 
patches. The lack of sharp spikes in the errors for the circular patches shows the effectiveness of applying angular 
dissipation. There is slightly more jaggedness near the equator for the circular vs square patch errors, but this is 
overbalanced by the relative smallness of the circular patch error. 



Variable 


circular patch 


circular, no dissipation 


square patch 


£L 2 (N) t=1 


2.05 


2.05 


2.05 


£l 2 (A0*=io 


2.05 


2.05 


2.04 


S L2 ( N )t=80 


2.04 


2.04 


2.01 


£l 2 (A0i=100 


2.01 


2.07 


2.05 


&„(JV) i=1 


2.04 


2.04 


2.04 


£ Loo {N)t =w 


2.04 


1.99 


2.04 


£ Lx ,(N) t = 90 


2.01 


2.01 


2.06 


£ Loo (N) t =ioo 


1.98 


2.00 


1.93 



TABLE VI: Convergence rates of the error £(N), measured at t = 1, t = 10, t — 90, and t = 100. 



Variable 


circular patch 


circular, no dissipation 


square patch 


£L 2 (iV*)t=i 


2.11 


2.10 


2.11 


£l 2 (iV* )t=io 


2.13 


2.13 


2.11 


£l 2 (Nw)t=90 


2.09 


2.09 


2.08 


£l 2 (Ny )t=ioo 


2.02 


1.98 


2.00 


£ ioc (iV*)t=i 


2.08 


2.08 


2.08 


&Loa (A r *)t=10 


2.09 


2.05 


2.10 


£l oc , (N\i)t=90 


2.05 


2.00 


2.06 


£l„ (A r *)t=ioo 


1.98 


2.01 


1.93 



TABLE VII: Convergence rates of the error £(Ny), measured at t = 1, t = 10, t = 90, and t — 100. 
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Variable 


circular patch 


circular, no dissipation 


square patch 




2.247 x 1(T 9 


3.325 x 1(T 9 


2.897 x 1CT 9 




1.706 x 10" 9 


2.747 x 10" 9 


2.315 x 1(T 9 



TABLE VIII: The values of the errors in JV and iV* measured at t = 100, 





FIG. 11: Surface plot of the error in N at t = 100 for the circular patch. 



FIG. 12: Surface plot of the error in N at t = 100 for the square patch. 




FIG 



13: 



Surface plot of the error in iV* at t — 100 for the circular patch. 



FIG. 14: Surface plot of the error in at t = 100 for the square patch. 
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VII. CONCLUSION 

We have proposed two new methods for enhancing the accuracy of CCE, One is a numerical method modifying the 
stereographic patches used in the characteristic evolution code to conform to the circular patch boundaries as used in 
meteorology [28|]. The other is a geometrical method that bases the waveform on the limiting behavior at T + of the 
Weyl tensor component ^4 rather than the news function. 

We have used a scalar wave testbed to compare the circular patch method against the cubed sphere method, which 
is also extensively used in meteorology (2t|. We found, for equivalent computational expense, that the cubed-sphere 
method has an edge in accuracy over the stereographic method. The cubed-sphere error in the scalar field £($) is » | 
the stereographic error but that the advantage is smaller for the higher derivatives required in gravitational waveform 
extraction. The cubed-sphere error £(95 2 <I>) is only s» | the stereographic error. An advantage of the stereographic 
approach is its relative programming simplicity. But as originally pointed out in [29j, and demonstrated recently for 
the case of a characteristic evolution code [33] , once all the necessary infrastructure for interpatch communication is in 
place, the shared boundaries of the cubed-sphere approach admit a highly scalable algorithm for parallel architectures. 

We used the circular patch stereographic code to compare waveform extraction in a linearized wave test directly via 
the Bondi news function N and its counterpart N-q, constructed from the Weyl curvature. For this purpose, we were 
able to successfully implement a new form of angular dissipation in the characteristic evolution code, which otherwise 
would be prone to high frequency error introduced by the irregular way a circular boundary cuts through a square 
grid. Our test results show that this dissipation works: the resulting error in the waveforms and metric quantities is 
smooth. In addition, the extensive analytic and numerical manipulations carried out to compute the limiting behavior 
of the Weyl curvature was demonstrated to yield second order accurate results for Ny. 

In the linearized tests presented here, neither N nor N-q, was a clear winner. We already knew that the original 
news module based upon a square stereographic patch worked well in the linear regime. The news N calculated on a 
circular patch had lower error than that on a square patch but only by a ~ 30% factor. In turn, the news calculated 
via on the circular patch had a lower error than N on the circular patch by a ss 24% factor. Weyl tensor extraction 
is slightly more accurate than news function extraction, even though there are many more terms involved. 

All errors were second order convergent. However, while there was a small fractional error « .1% in metric quantities 
such as J, the corresponding averaged error in the and N was « 4% for the circular patch runs and the maximum 
error at the equator was « 9%. These errors did not vary appreciably (w 30%) with the choice of discretization 
method, i.e. circular patch or square patch. They reflect the intrinsic difficulty in extracting waveforms due to the 
delicate cancellation of leading order terms in the underlying metric and connection when computing 0(1 /r) quantities 
such as $4. The excellent accuracy that we find for the metric suggests that perturbative waveform extraction must 
suffer the same difficulty. In that case it is just less obvious how to quantify the errors. The delicate issues involved 
at 2 + have been shown to have counterparts in extraction on a finite worldtube [24j . 

Waveforms are not easy to extract accurately. However, the convergence of our error measurements is a positive sign 
that higher order finite difference approximations might supply the accuracy that is needed for realistic astrophysical 
applications. Whether the advantages the new methods proposed here prove to be significant will depend upon the 
results of future application in the nonlinear regime. 
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